module mod_chem_species

   use iso_fortran_env, only: dp => real64
   use mod_chem_thermo
   use mod_chem_transport
   use stdlib_string_type

   private
   public :: type_chem_specy
#include <chem.f90>
   type :: type_chem_specy
      character(len=:), allocatable :: name
      character(len=:), allocatable :: geometry
      integer, allocatable :: composition(:)
      real(dp) :: molar_mass
      type(type_chem_thermo) :: thermo
      type(type_chem_transport) :: transport

   contains

      procedure :: initial => initial_species

      procedure :: get_cp_R => SpeciesGetCPR
      procedure :: get_cp_mole => SpeciesGetCPMole
      procedure :: get_cp_mass => SpeciesGetCPMass
      procedure :: get_enthalpy_RT => SpeciesGetEnthalpy_RT
      procedure :: get_enthalpy_mole => SpeciesGetEnthalpyMole
      procedure :: get_enthalpy_mass => SpeciesGetEnthalpyMass
      procedure :: get_internal_energy_RT => SpeciesGetIntEnergy_RT
      procedure :: get_internal_energy_mole => SpeciesGetIntEnergyMole
      procedure :: get_internal_energy_mass => SpeciesGetIntEnergyMass
      procedure :: get_entropy_R => SpeciesGetEntropyR
      procedure :: get_gibbs_RT => SpeciesGetGibbsRT
      procedure :: get_viscosity => SpeciesGetViscosity
      procedure :: get_thermal_conductivity => SpeciesGetThermalConductivity

   end type type_chem_specy

contains

   subroutine initial_species(this)
      class(type_chem_specy), intent(inout) :: this

      call this%thermo%initial(this%geometry)

   end subroutine initial_species

   function SpeciesGetCPR(this, T) result(cp)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: cp

      cp = this%thermo%get_cp_R(T)

   end function SpeciesGetCPR

   function SpeciesGetCPMole(this, T) result(cp)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: cp

      cp = this%get_cp_R(T)*R0

   end function SpeciesGetCPMole

   function SpeciesGetCPMass(this, T) result(cp)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: cp
      cp = this%get_cp_mole(T)/this%molar_mass
   end function SpeciesGetCPMass

   function SpeciesGetEnthalpy_RT(this, T) result(H_RT)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: H_RT

      H_RT = this%thermo%get_enthalpy_RT(T)

   end function SpeciesGetEnthalpy_RT

   function SpeciesGetEnthalpyMole(this, T) result(h)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: h

      if( T >=298.15_dp) then
         h = this%get_enthalpy_RT(T)*R0*T + this%thermo%h0
      else
         h = this%get_cp_mole(T) * (T-This%thermo%T0) &
             + this%get_enthalpy_RT(this%thermo%T0)*R0*T + this%thermo%h0
      endif

   end function SpeciesGetEnthalpyMole

   function SpeciesGetEnthalpyMass(this, T) result(h)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: h
      h = this%get_enthalpy_mole(T)/this%molar_mass
   end function SpeciesGetEnthalpyMass

   function SpeciesGetIntEnergy_RT(this, T) result(u_rt)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: u_rt
      u_rt = this%thermo%get_enthalpy_RT(T) - 1.0_dp

   end function SpeciesGetIntEnergy_RT

   function SpeciesGetIntEnergyMole(this, T) result(u)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: u

      u = this%get_enthalpy_mole(T) - R0*T

   end function SpeciesGetIntEnergyMole

   function SpeciesGetIntEnergyMass(this, T) result(u)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: u

      u = this%get_internal_energy_mole(T)/this%molar_mass

   end function SpeciesGetIntEnergyMass

   function SpeciesGetEntropyR(this, T) result(S_R)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: S_R
      S_R = this%thermo%get_entropy_R(T)
   end function SpeciesGetEntropyR

   elemental function SpeciesGetGibbsRT(this, T) result(G_RT)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: G_RT
      G_RT = this%thermo%get_gibbs_RT(T)
   end function SpeciesGetGibbsRT

   elemental function SpeciesGetViscosity(this, T) result(mu)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: mu
      mu = this%transport%get_viscosity(T)
   end function SpeciesGetViscosity

   elemental function SpeciesGetThermalConductivity(this, T) result(k)
      class(type_chem_specy), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: k
      k = this%transport%get_thermal_conductivity(T)
   end function SpeciesGetThermalConductivity

end module mod_chem_species
